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The dynamical mean field theory (DMFT), which is successful in the study of strongly correlated fermions, 
was recently extended to boson systems [Phys. Rev. B 77 , 235106 (2008)]. In this paper, we employ the 
bosonic DMFT to study the Bose-Hubbard model which describes on-site interacting bosons in a lattice. Using 
exact diagonalization as the impurity solver, we get the DMFT solutions for the Green's function, the occupation 
density, as well as the condensate fraction on a Bethe lattice. Various phases are identified: the Mott insulator, the 
Bose-Einstein condensed (BEC) phase, and the normal phase. At finite temperatures, we obtain the crossover 
between the Mott-like regime and the normal phase, as well as the BEC-to-normal phase transition. Phase 
diagrams on the fi/U — t/U plane and on the T/U — i/U plane are produced (f is the scaled hopping amplitude). 
We compare our results with the previous ones, and discuss the implication of these results to experiments. 

PACS numbers: 71.10.Fd, 67.85.Hj, 03.75.Hh, 05.30.Jp 



I. INTRODUCTION 

The ultracold atoms trapped in an optical lattice have 
aroused growing interests in recent years. By regulating the 
various parameters of the standing wave laser fields that create 
the optical potentials, such as the laser power and wave length, 
many theoretical models in the condensed matter physics can 
be realized experimentally, especially those for the strongly 
correlated many body systems. 1 In particular, bosons in a lat- 
tice have been widely studied in theory and experiment. The 
investigation of the correlated bosons can be traced back to 
the study of 4 He^ Recently D. Jaksch has pointed out that 
the Bose-Hubbard model (BHM) Eq.([T]i can well describe the 
ultracold boson atoms in an optical lattice, if one assumes a 
short-range pseudo potential interaction between the atoms 
and that the Wannier functions are well localized on the lattice 
site^i 

The BHM has been studied using various analytical and nu- 
merical methods. In their seminal paper, M. P. A. Fisher et 
al. used field theoretical approach to investigate the ground 
state of this model and obtained the superfluid (SF)-Mott in- 
sulator (MI) transition on the mean field levels The Mott in- 
sulator is an incompressible state where integer number of 
bosons are localized on each site, while the superfluid phase 
is compressible and has nonlocal boson wave functions. As 
an interesting phase, the ground state of MI is considered 
as a good candidate to realize the qubits for quantum in- 
formation processings 6 ^ For weakly interacting bosons for 
which the fluctuations around the mean-field state are small, 
the Bogolubov theory or the Gross-Pitaevskii equation applies 
but both fail to predict the SF-MI transition. 8,9 Beyond the 
mean field level, methods that can tackle the strong correla- 
tions have been applie d 3,5,9,10,11 ' 12 , including the Gutzwiller 
approac h 13,14 , Bethe ansatzi 5 -, time-dependent variational 
principle method 1 —, slave boson approac h 17 ' 18,19 , the strong 
coupling expansio n 20 ' 2 1 ,22,23 , variational method based on 
mean field theory 2 !, and the effective action approac h 24 ' 25 , 
etc. Recently the cavity method based on the Bethe lattice is 
also applied to BHM 2 ^. The numerical tools such as quantum 
Monte Carlo 27 ' 28 ' 29 ' 30,31 ' 32 ' 33 and density matrix renormaliza- 
tion group 34 ' 35 are used frequently for the unbiased studies. 



The physics of the BHM depends critically on its spatial di- 
mension. For a one dimensional system, the Mermin- Wagner 
theorem 36 excludes the possibility of Bose-Einstein Conden- 
sation (BEC) at finite temperatures. In the strong interaction 
regime, the bosons behave as the Tonks-Girardeau gas whose 
properties are similar to the noninteracting fermions i 4 ' 37 ' 38 
The BEC in two dimensions can be viewed as a quasi- 
condensate, and the Mott transition is of the Kosterlitz- 
Thouless type^ The MI is strictly defined at zero tempera- 
ture and the MI-BEC transition is a quantum transition well 
defined at zero temperature. However, there is a finite tem- 
perature range where the occupation is fixed at an integer. As 
temperature increases, the system changes from this Mott-like 
regime into the normal phase through a smooth crossover. For 
D > 2, the SF phase also transits into the normal phase at a 
finite critical temperaturei^iiS 

Experimentally, the BHM has been realized in the system 
of alkali metal atoms in an optical lattice. The SF to MI tran- 
sition has been observed in ID, 2D and 3D, by changing the 
depth of the optical lattices i 40 ' 41 ' 42 These studies focused on 
the MI and the SF phase mainly in 2D and 3D 39 ' 41 ' 43 ' 44 ' 45 at 
finite temperatures. The temperature in these studies is a key 
factor and its effects on the observation cannot be ignored. Re- 
cent studies show that due to the finite temperature effects, the 
sharp peaks in the momentum distribution commonly adopted 
to identify the SF cannot be used as a reliable signature! 27 ' 46 ' 47 
Due to technical difficulties, accurate studies for high dimen- 
sional system (D > 2) at finite temperatures that match the 
experiments are still highly desirable. 

Here, we are interested in the strongly correlated bosons 
on a high dimensional lattice (D > 2) at finite tempera- 
tures, which received less attention in previous theoretical 
studies. One suitable method for our purpose is the dynam- 
ical mean field theory (DMFT). DMFT is an exact theory in 
infinite spatial dimensions. 48 As an approximation for finite 
dimensional systems, it has been widely used in the study 
of strongly correlated fermion systems and received much 
success! 9 - In a recent work, K. Byczuk and D. Vollhardt ex- 
tended the idea of DMFT to correlated boson systems and ap- 
plied it to the bosonic Falicov-Kimball model.— In their the- 
ory, instead of the usual way of scaling the hopping amplitude 
in the fermionic DMFT, they used a different scaling ansatz 
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for bosons: scaling f,y — > fy/z' 1- - 7 ' for terms containing the 
anomalous averages (£>) or {b'~}, and scaling — > f, ; 7 -y5' ! ~ ; ' 
for others. Such a new scaling is used to guarantee that the 
energy density is finite as the spatial dimension goes to in- 
finity, even if anomalous averages are involved in the boson 
systems. It should be noted that the derivation of the above 
B-DMFT equations is not unique. Recently, the essentially 
identical equations are also obtained by using an uniform 
scaling of t — > t/z and keeping up to the subleading order 
in the 1/z expansion. 51 In this paper, we apply B-DMFT to 
the BHM. The resulting effective bosonic impurity model is 
solved by exact diagonalization (ED) method 50 ' 52 . Results are 
presented for various phases at finite temperatures and com- 
pared to other theories and the experiments. 

This paper is organized as follows. In Sec. II, we briefly in- 
troduce the single band Bose Hubbard model. We present the 
B-DMFT equations for the BHM, detail the impurity solver 
that we use, and give tests and benchmarks. In Sec. Ill, the 
main results of B-DMFT are shown and discussed. Sec. IV is 
a conclusion. We put some technical details in Appendices. 



II. MODEL AND METHOD 
A. Bose Hubbard Model 

The single band BHM is defined by the Hamiltonian below 

(i.j) i » 

where b\ and b, are the boson creation and annihilation op- 
erators on site i, respectively. They obey the commutation 
relation bj] — 6,j. «, = b\bi is the boson number operator 
on the site i. Here, we consider the hopping amplitude = t 
which is nonzero only for the nearest neighbors and U is the 
on-site energy. Feshbach resonances can be used to change 
the interaction strength over a wide range, even from repul- 
sive to attractive. 4 In this paper we study the repulsive BHM 
with U > 0. We add the chemical potential which controls 
the number of bosons in the grand canonical ensemble. For 
simplicity, we take the density of states of the Bethe lattice, 
which is semicircular in the limit of infinite coordinations, 

D(e) = -L V4f 2 - e\ (|e| < It). (2) 



B. B-DMFT Equations 

In DMFT, a lattice model is mapped into a single impu- 
rity problem with the self-consistently determined bath spec- 
tra. It becomes exact when the spatial dimension is infinite 
and hence ignores the spatial fluctuations from the outset. 
However, it fully takes into account the temporal fluctuations 
(imaginary time)42i The key ingredient to extend the DMFT to 
bosons is a proper scaling of the hopping amplitude of bosons 
in the limit of infinite dimensions^ 



We adopt the ansatz of scaling in Ref. and imple- 
ment similar derivations for the BHM. The detail of deriva- 
tions can be found in Ref. HO- Here we present only the 
final B-DMFT equations. For simplicity, here we use the 
Nambu representation^ for the boson operators b + (r) 
j b^(r) b(r) ), and for the on-site interacting Green's func- 
tions (GFs) as in Ref. 

G(r-r') = -<r T [b(T)b t (T')]> 

= / -{T T [b{r)b\r')]) -{T T [b{T)b{r')]) \ 
\ -{T T [bHT)bHr')]) -(T T [b\ T )b(T')]) J (3) 

_(G 1 (t-t') G 2 (t-t')\ 
\G 3 (t-t') G 4 (t-t') J' 

According to the definition, the following relations hold for 
the components Gt,(t - r') = G* 2 (t - t') and G^{t - t') = 
q(r-r'). 

The action for the effective impurity model obtained 
through the cavity method 4 - reads, 

S *ff = if dT if dr'bJ(T)[-^- 1 (T-T')]bo(T') 

+ Jf drf n (T)[no(r) - 1] + jf drOjb (r). (4) 

In this equation, Oo is related to the condensation via 

O; = ( -t{bl)s,„ -t(b ) Seff ). (5) 

Here t is the hopping amplitude after the scaling has been car- 
ried out. (b^) is treated as a r-independent quantity since we 
are studying an equilibrium theory. The Weiss field Q^. (icj n ) 
represents the effective field from the environmental fluctua- 
tions acting on the impurity site. 

The self energy is defined through the Dyson equation 

S(/w„) = 20Q l (io) n ) - G~\ia)„). (6) 

Here, G c is the connected GF defined as 

G c (t - r') = G(t - t') - G dis (r - t'), (7) 

where Grf, s (r - t') is the disconnected part. Its fourier trans- 
form is given in Appendix C. In the imaginary time axis it is 
a constant and coincides with the condensed fraction in the 
thermal dynamical limit. £ and Qq have the same symmetry 
properties as the GF. Among the four matrix elements only 
two functions are independent. It is noted that the defini- 
tion of the self-energy in Eq.© has a factor of 2 difference 
from Eq.(l 1) in Ref. [50i This difference can be traced back 
to the different conventions used for path integrals in Nambu 
representation! 5 ^ We have checked that our equations are self- 
consistent and they guarantee £(/«„) = for U = 0. 

The connected local GF of the lattice Hamiltonian is given 
by the lattice Dyson equation 

G c (ioj n ) = ^-Y [[Gf ]-'(*, io> n ) - nico n )\ X . (8) 
Mum *y 
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Here Ni alt is the total lattice number. Gf (k, ia>„) is the con- and 



nected GF of the non-interacting system Ho = 
It reads 

Gf\k, iu)„) = [ia>„o- 3 - (e k - fi)1]~ 



In the actual calculations, we transform the summation over 
k in Eq.® into the integral over energy. The explicit integral 
formulas involving the semicircular density of states Eq. (lA2b 
are summarized in Appendix A. Eq.(f3|-Eq.((9| constitute the 
B-DMFT self-consistency equations for the BHM. 



C. Impurity Solver 

In order to solve the B-DMFT equations, a suitable impu- 
rity solver should be selected. To avoid technical complexities 
we use the exact diagonalization method to solve the impurity 
model. It is simple, fast, while at the same time qualitatively 
keeps the nontrivial many -body physics of the problem^. The 
effective impurity Hamiltonian equivalent to the action Eq.© 
reads 



k=l k=\ 

U t 
+—n (no - 1) + %b . 



(10) 



The creation and annihilation operators a! and & k are for the 
environmental degrees of freedom and are all in the Nambu 
representation. B s is the number of bath sites. Ea and V* 
are the kinetic energy of environmental bosons and the cou- 
pling strength between the environment and the impurity, re- 

E k \ E k 2 



spectively. They are 2x2 matrices, Ea = 

'V a V k2 



Vt = 



V a V k4 



Eki E k 4 

From the Hermiticity of H mp - H 



and 

imp'' 

and 



we have E k4 = E k \ being real, E k3 = £* 2 , V k3 = V* k2 
V k 4 = VTj . The requirement that Hj mp is equivalent to the ef- 
fective action S e /f in Eq.® gives the following relation (see 
Appendix B) between £? 1 and E*, Va, 



-lCO n (T 3 + 



(11) 



We solve the B-DMFT equations using an iterative scheme 
as usually done for fermions. We start from an initialization 
of the parameters Ej, Y k (k — 1, ..,B S ) and Oq. With them 
we calculate @ l (iu>„) and define the impurity model Eq.dTOb. 
The impurity Hamiltonian is then solved by ED to produce the 
connected GF G e and a new Oo, according to the following 
equation, 



4»o 



-?<bo>. 



(13) 



(9) Here (...) represents the average under Hj mp . G is calculated 
from the Lehmann representation. Grf, s (/w„) is the discon- 
nected Green's function. Details are in Appendix C. 

Using Eq.©, one obtains the self-energy "L from G c . It is 
then put into the lattice Dyson equation Eq.© to produce a 
new G c . Using Eq.© again, we update the Weiss field @q 
and from it we get the new parameters E k and Va through the 
following fitting procedure.^ A distance function D [E k , Va] 
is defined as 



£>[Ea,Va] 



2 



\@0l(iten) 



@ 2(i(O n )\ 



(14) 



Here @ Q l is calculated from Ea and Va through Eq.(fTT1). 
D [Ea, Va] is then minimized with respect to Ea and Va to find 
the optimal parameters that can reproduce 0^ x {i(jJ„) best. With 
these optimal parameters we define a new impurity model to 
be diagonalized again. The iteration continues until the lattice 
GF converges. 

One specialty of boson is that it has an infinitely large local 
Hilbert space. This poses difficulty for ED-based numerical 
methods when adapted for bosons^ In our ED calculations, 
we truncate the local Hilbert space by using 1 boson states 
for each mode, with a finite number. As the simplest algo- 
rithm, we use the boson number eigen state \n) (n = 0, 1, N) 
as our local basis, and keep the implementation of optimal 
basis for a future improvement M The truncation of boson 
Hilbert space introduces additional approximation and influ- 
ences the accuracy of our results, especially in the BEC phase 
(see below). It is therefore important to check our results with 
respect to and to make sure that the truncation errors are 
under control. 

However, the truncation described above introduces a new 
problem to the commutation relation of boson operators. In 
the truncated Hilbert space, one has 
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V2 
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o V^v 





(15) 



G c (/w„) = G(ia) n ) - G dis (icj n ), 



(12) 



and £>' is the hermitian conjugate matrix of b. From 
these one gets W?' = diag{\,2, ...,N, 0) and b^b = 
diag{0,l, ...,N - l,N}. The commutation relation reads 
[b, b^] = diag{ 1,1,..., -AO with the incorrect trace Tr[b, b^] = 
0. Therefore, using representation Eq.dl3t in our calculation 
will lead to incorrect weight in the GFs as well as in the den- 
sity of states. This problem cannot be remedied by increasing 
N. To overcome this difficulty, we modify the representation 
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of b (b 1 accordingly) into 



b = 
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Vjv+ i J 



(16) 



It produces bb = diag{0, 1, ...,N- l,2N + 1}, and [bb']a = i 
(i = 1, ...,#+ 1) and [bb%, N+1 = [bb% +hN = y/N(N + 1). 
The trace Tr[b, b] — is still incorrect. However, b'b from 
representation Eq.dT3]> and bb from Eq.dToTi, if combined to- 
gether, give the correct trace Tr[b, b] = N + 1 Therefore, 
our strategy is that for any operators involving bb, such as 
(i\b T \j)(j\b\i) in the Lehmann representation of the diagonal 
GF, we use Eq.(fT31). For operators involving bb' such as 
{i\b\j){j\tf\i), we use Eq.® (see Eq.(C2) and (C3) in Ap- 
pendix C). In this way, the truncation introduced boson com- 
mutation problem is solved. 

For the bosonic impurity model Eq.([T0l) with B s bath sites 
and N + I states for each boson mode, the dimension of the 
Hilbert space is S , — (N + l) Bs . To describe the BEC phase 
where Oo + 0, the total particle number can no longer be 
used as a good quantum number. In this case both ED and 
calculating the GFs are very time consuming. As a result, 
the parameters N and B s are severely limited by the present 
computer power. In the DMFT (ED) study of the Fermi Hub- 
bard model with four states on each site, it was shown that 
results converge quite fast with the bath site number B s , and 
B s = 4 already gives qualitatively reliable results^ To ex- 
plore the B s and N dependence of calculations for the Bose 
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FIG. 1: (Color online) The total boson occupation as functions of 
chemical potential yu. (a) and (b): U = 50.0, t = 0, and T = 1.0; (c) 
and (d): U = 0, t = 1.0, and T = 1.0. 



Hubbard model, we calculate the N, ol - ft curves in both the 
atomic and the free boson limit for different and B s values. 
The results are shown in FigQ] In both limits, as long as the 
bath site number B s is larger than zero, the results are already 
very close to the exact ones. At the same time, A^-dependence 
is more severe. The curves keep improving observably until 
N >\Q. Taking a compromise between Af and B s , in our study 
we do all the calculations at B s = 1 (one bath site) and N = 15 
(16 boson states) unless stated otherwise. We check our re- 
sults using larger Af and B s and make sure that our conclusion 
does not depend on the selection of B s and N. For the Mat- 
subara frequencies co„ = Inn/B in the GFs, we use the cut off 
\u„\ < 2000. 



D. Non-Interacting Limit and Atomic Limit 

In this section, we check the B-DMFT formulas and our 
numerical results in the noninteracting as well as in the atomic 
limit. For a noninteracting boson system, the bosons can move 
freely in the lattices. They condense into a single particle state 
when the temperature is lower than T c . The exact solution 
of the Bose Hubbard model in this limit gives the thermally 
excited boson occupation N e as 



de- 



Die) 



P(e-M) - 1 ' 



(17) 



The B-DMFT equations Eq.(f3)-Eq.((9]l can also be solved ex- 
actly in this limit. 50 At U = 0, by carrying out the Gaussian 
integral in Eq.© and doing the functional derivative of the 
free energy with respect to £? ' (subtracting the disconnected 
contribution), we get the connected GF as 



G c (ia) n ) = -Qo(ioJ„) 



(18) 



which is independent of <E>o. This gives a zero self-energy 
T,(ici>„) = according to Eq.©, as expected. From Eq.© and 
(O we obtain the GF 



AeD(e) [ico„(T 3 - (e - yu)I]~ 

OO 



(19) 



It is the exact result for the free bosons. The expression 
Eq.(TT7l for the thermal excited particle number can be re- 
covered from it using the fluctuation-dissipation theorem. 
The order parameter of BEC reads 



= TOjgo(iO). 



Together with Eq.©, it gives 



4>o 



= 0, fi < -It, for normal phase, 
* 0, yu = -It, for BEC phase. 



(20) 



(21) 
(210 



Considering that for fj. = -2t, G c (/0) = -(l/f)I, we cannot 
determine the nonzero value of Oq in the BEC phase solely 
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from the self-consistency equations. This corresponds to the 
fact that for free bosons, the condensed fraction cannot be 
fixed without giving the total particle number N tot . The re- 
sults above can also be obtained from the equation of motion 
of GFs, starting from the impurity Hamiltonian Eq.dTOb with 
U = 0. From this approach one gets 



G c (ia>„) = 



iu>„<T 3 + /il ~ 4 V V* (iaj„cr 3 - 2E k ) 1 V] 



k 



(22) 



number of localized particles at zero temperature. As the tem- 
perature increases, thermal fluctuations dominate and the lo- 
calized state will melt gradually. In this limit, <J>o = ac- 
cording to Eq.©. The density of states Z?(e) becomes a delta 
function, and Eq.® reduces to 



G c (ia>„) = [m n cr3 T.(ia) n )]~ 



(24) 



Comparing with Eq.©, one gets @ l (iu>„) = (icj n cri, + yul) /2. 
When inserted into the effective action Eq.©, it gives ex- 



and 



(bo 1 ') = ®J 



fi + 2y,v i E: l v] 



k 



(23) 



Substituting the parameters and V* with the Qo in Eq.dTTI). 
we can get the same results as Eq.dTSb and d20b . 

In FigfJl the real and imaginary part of the self-energy are 
shown for U — 0. Both the diagonal and off-diagonal compo- 
nents tend to zero as the number of boson states N increases. 
The self energy is not strictly zero for finite N, because in the 
truncated space, boson operators do not obey canonical com- 
mutation relations and even U = does not correspond to a 
free system. From Fig(3j it is seen that the N e - fi curves from 
B-DMFT at U = agree well with the exact ones. The devia- 
tion at high temperatures decreases as N increases, consistent 
with what we find in Fig [2] 

In the atomic limit t — 0, the quantum fluctuations of the 
boson number operators disappear. Each site has an integer 
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FIG. 3: (Color online) DMFT result (dots with eye-guiding lines) 
and the exact result (solid line) for the thermal excited boson number 
N e at V = 0, t = 1.0. (a) as a function of /j at T = 1.0; (b) as a 
function of T at fi = —21. Insert: N e as a function of N at /j = —2t, 
T = 1.5. 




FIG. 2: (Color online) The self-energy of the non-interacting system 
as a function of Matsubara frequency, (a) and (b): diagonal compo- 
nent Si; (c) and (d): off-diagonal component £2- They are calculated 
at U = 0, t = 1,0, ju = —2.1, and T = 1.0. Symbols are denoted in 
the figure. 



FIG. 4: (Color online) The total number of bosons as functions of 
the chemical potential /i in the atomic limit with U = 50.0 for dif- 
ferent temperatures. The lines are the exact results and symbols for 
B-DMFT results. Inset: the compressibility dN tot /dfl as functions of 
yL obtained from B-DMFT. 
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actly the action in the atomic limit. Our numerical results 
for the atomic limit obtained using B s — 1 and N - 15 are 
shown in FigH] The B-DMFT results (squares with guiding 
lines) agree well with the exact ones (curves). Note that in 
the atomic limit, the B-DMFT results depend very weakly on 
B s and N. It is seen that the thermal activation will smear the 
Mott plateaus and the compressibility dN lo ,/d/-i has broadened 
peaks. For U = 50, the Mott plateaus are clear at low tempera- 
tures and their features disappear completely at about T = 10. 
This observation agrees with the conclusion that the MI melts 
completely at about T* = 0.2U in the limit t = 0^ 



in. RESULTS AND DISCUSSIONS 

In this section we discuss the physical results obtained by 
the B-DMFT for the BHM. At zero temperature, the system 
should be either in the BEC phase for weak interaction or in 
the MI phase for strong interaction. At finite temperatures, 
besides the BEC and MI phases that are extended from the 
ground state, there is the normal phase that is connected to 
the MI and BEC phase in low temperature regimes, through a 
crossover and a second-order phase transition, respectively. 
FigEJa) and (b) show the diagonal component of the con- 
verged connected GFs typical for the BEC, MI, and normal 
phases. They are calculated at a finite but low temperature 
T/U = 0.05. All the high energy parts show the expected 
behavior ReG f i(/a>„) °c l/o% and ImG e i(/ai„) °c 1/oj„. The 
low energy behaviors are markedly different between the MI 
phase and the other two. Similar to the MI phase of fermions, 
ImG c i(za> n ) — > at zero frequency, signaling strong scatter- 
ing and nonexistence of well defined low energy quasiparti- 
cles in the MI phase. In contrast, in the BEC and the normal 
phases, the diagonal connected GFs are qualitatively similar. 



In Figj5jc) and (d) are the corresponding off diagonal GFs in 
the three phases. The condensation in the BEC phase con- 
tributes to a sharp peak in the low energy regime of ReG C 2. 
While in the MI and the normal phases, no such peak appears. 
In all the three phases, ImG C 2 — as it is required by its defi- 
nition and symmetry. 

The full GF (not shown in PigO can be written as the form 
G m (i(jL>„) = G cm (iu> n ) + A m 6 n fi, (m = 1,2). For the BEC phase 
shown in Fig|5] our numerical calculation gives Ai = -2.72, 
A2 = -2.71, and -B(bo) 2 = -2.61. Within numerical errors, 
our result is consistent with the equation Aj = A2 = -B(bo) 2 
as can be seen from the Lehmann expressions for GF in Ap- 
pdix C. For the MI and the normal phases, we always get 
Ai = A 2 = 0. 

The total particle occupation is calculated by 

N tot = ~\Yj G d(i"n)e ibJ * * + (b ) 2 , (25) 

where the condensed boson number A^o is 

A'o = <fto> 2 . (26) 

For U > 0, a simple mean field analysis shows that the free 
energy becomes a quartic function for large (bo) and hence 
both A^o and N tot can be determined solely by the B-DMFT 
equations. They are plotted in Fig|6]as functions of the chem- 
ical potential. In Fig|6|a) we fix T/U = 0.05 and study the 
evolution of the curves as t/U decreases. For large t/U (small 
U for fixed t), N tot and A^o are increasing functions of /j, up 
to a boundary of fi and the system always stay in the BEC 
phase. At the boundary, the convergence becomes slow and 
difficult. As t/U is smaller, a plateau of N tot = 1 begins to 
emerge in the Ntot ~ f 1 curve. A'o has a temporal decreases 
at the corresponding fi and then continues to increase. The 
system is still in BEC phase, but the plateau and the dip in A^o 




-30 -20 -10 10 20 30 -30 -20 -10 10 20 30 

FIG. 5: (Color online) The connected diagonal GF G c \ ((a) and (b)) 
and off diagonal GF G C 2 ((c) and (d)) in the three phases: BEC phase 
(circle): t/U = 0.2; normal phase (triangle): t/U = 0.11; MI phase 
(pentacle): t/U = 0.05. All are calculated at t= 1.0, fi/U = 0.5, and 
T/U = 0.05. 
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FIG. 6: (Color online) N, ot and A'o as functions of /j. (a) and (b) 
T/U = 0.05 for different t/U (t = 1.0); (c) t = 1.0, t/U = 0.05 for 
different temperatures. Symbols with eye-guiding lines are denoted 
in the figure. 
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show the precursor to the MI phase. As tj U still decreases, the 
plateau enlarges at N m = 1 and the next one at N tot = 2 begins 
to appear, forming Mott-like regimes. No has a well formed 
gap corresponding to each plateau, and has a peak signaling 
BEC between two neighboring gaps. These BEC phases ap- 
pear around fi = 0, U, 2U, ... where two adjacent Mott plateaus 
are connected. As t / U decreases, the height of the A^o peak 
decreases and the critical temperature T c also decreases (see 
Fig|8jb)). For very large U such as t/U = 0.02, BEC doesn't 
appear any more because the critical temperature is lower than 
the actual T, (T/U) c < T/U = 0.05 at t/U = 0.02. FigHc) 
shows the temperature evolution at a fixed t/U = 0.05. The 
temperature effects on the Mott plateau as well as on the BEC 
phase are clearly observable. As temperature rises, the Mott 
plateaus gradually blur at the shoulders and the condensed bo- 
son number A^o reduces to zero. For high enough temperature, 
the Mott plateaus finally disappear and the MI crosses over to 
the normal phase. 

A recent B-DMFT study for the bosonic Falicov-Kimball 
model reveals that the local repulsion enhances the transition 
temperature of BEC. 5l i Here we study the influence of U on 
BEC in the BHM. In Fig0 we plot the N - T curves at dif- 
ferent t/U values for a fixed total density N, , = 1.5. For a 
given t / U, No is a decreasing function of T, and reduces to 
zero at T c . With decreasing t/U values (increasing U/t), the 
No - T curve shifts downwards, leading to smaller A^o for a 
given T and a reduction of T c . This is consistent with the 
naive picture that strong local repulsion between bosons tends 
to suppress the particle number fluctuations and act against 
the BEC. Also, the quasiparticle states into which the bosons 
can condense are turned into incoherent states and shifted into 




FIG. 7: (Color online) The condensed boson number No as functions 
of temperature T for fixed total number N M = 1.5 and different t/U 
(f = 1.0). From top to bottom, t/U = 0.5,0.3,0.15,0.1,0.02, respec- 
tively. Inset: Nq changes with t/U for T = 0.3 (circle) and T = 2.0 
(pentacle), respectively. The dashed line is No = 1.5. Solid lines are 
guiding lines. 



the Hubbard bands by a large U. Therefore, our conclusion 
is that, different from the bosonic Falicov-Kimball model, the 
local repulsion in the BHM reduces the transition temperature 
of BEC. 

Different phases in the system can be distinguished from 
N tot and A^o- At zero temperature, the system has two phases: 
the BEC phase and the MI. Due to the competition between 
the on-site repulsion U and the hopping f, there is a quan- 
tum phase transition between them. As temperature increases, 
the BEC phase and the MI will change into normal phase 
through a phase transition and a crossover, respectively. We 
have therefore three phases to identify at finite temperatures: 
a nonzero No signals the BEC phase, while the MI has No = 
and an integer N ro , with zero compressibility dN, ol /dfi; the 
phase with No - but a finite compressibility is the normal 
phase. According to this criterion, we plot the phase diagrams 
inFigl 

In FigJSJa) is the phase diagram on the fi/U - t/U plane 
for finite temperature T/U = 0.05. It is obtained by scanning 
fi/U at fixed t/U. Three regimes are clearly shown, the MI 
phase, BEC and the normal phase. Due to the small N param- 
eter that we use, only the N tot = 1 and part of the N, ot = 2 MI 
domains are obtained. In the large t/U regime, BEC is stable. 
Between the two boundaries (circles and pentacles) is the nor- 
mal phase. The melting temperature T* and BEC transition 
temperature T c are marked by pentacles and circles, respec- 
tively. Similar finite temperature phase diagram is also obtain- 
able from a static mean field theory-^. To understand the tem- 
perature effects on this diagram, we resort to phase diagrams 
on the T/U - t/U plane at two different fi/U values, Fig|8jb) 
and (c). At T = 0, a BEC-MI quantum phase transition oc- 
curs at a critical t/U. The difference between FigJSJb) and (c) 
shows that (t/U) c is dependent on fi/U, consistent with the 




t/u t/u 

FIG. 8: (Color online) Phase diagrams, (a) in fi/U - t/U plane at 
T/U = 0.05; (b) and (c) in T/U - t/U plane at fi/U = 0.5 and 
fi/U = 1.0, respectively. Three phases, BEC, normal phase, and the 
MI phase are marked out in the figures. Lines are for eye-guiding. 



8 



lobe shape of the MI boundary in FigJHta). At finite tempera- 
tures, the normal phase appears as a quantum critical regime 
extending from the T — quantum critical point. This hints 
that in the normal phase near (t/U) c , critical behavior such 
as power law correlation should exist. Experimental observa- 
tion of such quantum critical features in the normal phase near 
(t/U) c will be an interesting issue. 

A representative quantity for comparison between differ- 
ent theories is the critical value (t/U) c at the tip of the the 
n — 1 Mott lobe. It has been obtained by various meth- 
ods. For the 3D cubic lattice, the world line QMC gives 
(t/U) c = 0.032(Ref. |H) and the worm algorithm QMC gives 
(t/U) c = 0.03408 (Ref. ©. Recent studies for the Bethe 
lattice with z = 6 give (t/U) c = 0.033 (Ref. H3) and (t/U) c = 
0.032 (Ref. ED). In our study, we obtained (t/U) c = 0.12 
for n/U = 0.5 (n=l Mott lobe). The apparent discrepancy 
between our value and the previous ones is because we didn't 
use the realistic lattice structure. In our calculations, we take 
z = oo literally and set t — 1 as the energy unit, hence z doesn't 
appear explicitly. Since different scalings relating t to t are 
used in B-DMFT for the normal and the condensed bosons, a 
critical value for z — 6 cannot be simply recovered from our 
result by doing an inverse scaling. A crude estimation, how- 
ever, gives (t/U) c ~ 0.12/z ~ 0.12/ Vz = 0.02 ~ 0.049 for 
z = 6, being consistent with the more accurate values. 

For realistic lattices with a finite coordinate z, the B-DMFT 
is still applicable but should be regarded as an approxima- 
tion to finite dimensional systems. For such a calculation, 
one should use the actual dispersion of the given lattice 
in Eq.©, and replace t with zt in Eq.©. Experimentally, 
the BEC-MI transition point was observed for &1 Rb atoms in 
3D optical lattices. 40 The transition occurs at a potential depth 
of l3E r which compares favorably with the mean-field value 
U/t = 5.8z j 9 ' 10 - 12 but differs from the more accurate QMC re- 
sults cited above. The B-DMFT calculation for the BHM in 
3D cubic lattice and quantitative comparison with the experi- 
ments as well as with the previous theoretical results will be 
an interesting topic. But this is only attainable when an ac- 
curate impurity solver is available. Therefore we leave it for 
future study. 

Finally we note that the ED method used in this work poses 
limitations to our study. Due to finite number of boson states 
N = 15, reliable calculations are only possible in the small 
N to , and small iVo regimes. The small number of bath sites 
B s = 1 causes slow convergence, especially near the MI- 
BEC transition. Therefore, for practical applications of B- 
DMFT to boson systems, it is necessary to develop an accurate 
and fast impurity solver. In this respect, the recently devel- 
oped bosonic NRG is a promising technique. 52 ' 56 For the ED 



method, an algorithm adopting the optimal boson basis will 
be interesting and progress is being made in this direction. 



IV. CONCLUSION 

In this paper we have performed the B-DMFT study for 
the BHM. Following the ansatz of scaling in Ref. we ob- 
tain the B-DMFT equations for the BHM. The bosonic ef- 
fective impurity Hamiltonian is solved by ED method with 
truncated boson Hilbert space. We focus on the finite tem- 
perature properties of the correlated bosons, and identify the 
MI, BEC and the normal phases. The repulsive U is found to 
suppress the BEC transition temperature T c . Phase diagrams 
on the fi/U - t/U and T/U - t/U planes are obtained, which 
disclose the quantum critical nature of the low temperature 
normal phase. Relevance of our results to other theoretical 
ones and the experimental observations are discussed. 
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APPENDIX A: INTEGRAL OF SEMICIRCULAR DENSITY 
OF STATES 

The lattice Dyson equation in the B-DMFT equations 
(Eq.©) is usually transformed into an integral over e of the 
form 

G c (i(L>„) = J deD(e)[uo n cr3 - (e-ju)I 

-2g \iu> n ) + [G»r\iun)Y l - (Al) 

One needs to carry out the integral for each ai n . For the semi- 
circular D(e) given in Eq.(2) 

D(e) = ^ V4f2 - e2 > (H * 2f ), (A2) 
the exact integral formula is given in the following. 
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de 



Die) 



' ^-Sgn(Im^) ^[e^Af- 
2t 2 

f-Sgn(Re^) ^fe^Afi 



ir- 



l 2f 2 ' 



In this equation, Im£ and Re£ are the imaginary and real parts 
of Sgn(x) = 1 and Sgn(x) = -1 for x being a positive and a 
negative real number, respectively. 



APPENDIX B: EFFECTIVE IMPURITY HAMILTONIAN 
AND ITS ACTION 

The statistical action for the impurity model Eq.dTOt reads 



Im£ * 0, 

Im£ = and |f| > It, 
Im£ = and |£| < It. 



(A3) 

(A3') 
(A3") 



Simp = dr ^ aJ(T) |ifl T tr 3 +E jfe ja*(T)+bJ(T)|i ( 9 T o-3 - ^/jb (r) 

+ dr V [aJ(T)V t bo(T) + bJ(r)V t a,(r)] + -„ (t)[« (t) - 1] + Oj(r)b (r) 

Jo [k=i L 

I 



(Bl) 



The partition function can be expressed as the path integral where Z a is the partition function of the bath degrees of free- 
over complex boson fields dom, and the effective action S b for the impurity is given by 



Z = f S>b;(T)S>b Q (T) f H $>al(T)2>a k (T)e- s "»". (B2) 

Carrying out the Gaussian integral for the environmental de- 
grees of freedom ak ' and ak, one obtains 



Z = Z 



(T)@b Q (T)e- s », 



(B3) 



' h = J u dr b5(T) -d T cr 3 _ -jj - g V* ^-3 T <r 3 + E*j V x 



bo(r) + |n (T)[no(r) - 1] + Oj(r)b (r) 



(B4) 



Comparing this equation with the effective action derived Through this equation the Weiss field Q^ 1 is related to the 

impurity parameters E* and Va. After a Fourier transform, 
one gets Eq. dTTb . 



from the cavity method Eq.©, one gets 
^o'(r-T') 

jd T a 3 - j/A - y V* |ia T cr 3 + E t ) V] 



ff(T - r'). 
(B5) 
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APPENDIX C: LEHMANN REPRESENTATION OF THE 
BOSON GREEN'S FUNCTION 



Here, Z = £]; exp(-f3Ej) is the partition function and ft — l/T 
the inverse temperature. The eigenvectors can be expanded by 
the basic vectors in the boson Fock space {|«)}, |i) = 2„ 



for w„ # : 



for u)„ — : 



The disconnected GFs Qdisii^n) are defined as the sum over 



They have the symmetric relation G^t — t') = G* 2 {t - r') and 
Gh(t — t') = G*(r - t'). Here we only consider G\ and G 2 . 



(C2) 



(C3) 



where (n\b\m) and (m|fc^|«) are the matrix elements discussed 
in Eq.(IT5b and Eq.dToTi. The coefficients A' n can be obtained 
from ED. 



(C5) 



(C6) 

I 

Ej = Ej parts in Eq.(|C3]) and dC6l . 



The boson GFs are calculated from their Lehmann repre- 
sentation after the eigenvalues and the eigenvectors are ob- 
tained by ED. In this appendix we present the corresponding 
Lehmann representations. The GFs are defined in Eq.© as 

G(T-T') = -<7V[b(T)bV)]> 

-<7V[fo(r)Z>t(T')]> -{T T [b(T)b( T ')]) \ 

-{T T [b\T)bHT')}) -{T T [b\T)b(T')}) J (CI) The diagonal GF is expressed as 

d(T-T') G 2 (r-r') 
G 3 (t-t") G 4 (t-t') 



J 



for oj„ + : 

1 x-i e~~P E i - e^ E ' 1 x-n e~ pE ' - e~P E > 

Gi(md = X . -7^ v A i\b\j){0\i) ~ X — — —Ai\b ] 

2Z t—t iai„ + (E, - Ej) 2Z itu„ + (Ej- En 

ij u 



\j)(j\b\i); 



for u>„ — : 

Gi(iO) 



^ Z F F (WM^IO ~ ^ E — = — = — <£[* + lj>0"I^IO 

2Z £^ £ >"~^ 2Z £^ ^~ £ < 

■A y e-^Kwx^io - 4 Z /;/ ./x//"->- 

I 



(i\b\j)=J]KAi(n\b\m), (C4a) 

(# + IO = Z A - A »< m l fot l n >- < C4b ) 

™ The off-diagonal GF reads 
I 

i „ e -PE, _ e -m 

Giiiun) = --Y- 7^ ^<WX#I0, 

Z ia>„ + (Ej — Ej) 
•j 

g 2 (/o) = ~2j <w)oifcio - 1 2j «*<wwi>- 
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